##
## ``The Rise Predicts the Fall"
##
## Jun Koga Sudduth and Curtis Bell
## Contact: jun.koga@strath.ac.uk
## April 2017
##
##


library(foreign)
library(MASS)
library(nnet)


##
## Data
##


dat <- read.dta("C:/Users/pjb13215/Dropbox/Entry & Exit manner in Dictatorships/FINAL SUBMISSION TO ISQ/ISQ_Sudduth_Bell_entryexit.dta" )

dat <- dat[order(dat$ccode, dat$leadercode, dat$year ),]


##
## Descriptive Statistics
##

names(dat)
dim(dat)  # [[1] 5814
length(table(dat$leadercode))  # 773 all leaders (including puppet)




##
## Entry manner

unique( dat[, c("entry"  ,   "entercode")])

#                           entry entercode
#            Successful Rebellion         1
#                       Protest           2
#                     Regime Coup         3
#           Foreign Installation          4
#                    Leader Coup         10
#   Irregular Election/Selection         11
#       First Election/Selection         12
#     Regular Election/Selection         20


entrydat <- dat[, c("ccode","stateabb"  , "leadercode", "leader" ,"entry"  ,   "entercode")]
head(entrydat)

head(unique(entrydat))

dim((unique(entrydat))) # 773 leaders

head(entrydat[entrydat$entercode==1,])
dim(unique(entrydat[entrydat$entercode==1,]))[1]    #[1] 27 successful rebellion entry leader
#unique(entrydat[entrydat$entercode==1,])   # check the list of rebellion entry leaders
dim(unique(entrydat[entrydat$entercode==2,]))[1]    # 4 protest entry leader
dim(unique(entrydat[entrydat$entercode==3,]))[1]    # 108 regime-coup entry leader
dim(unique(entrydat[entrydat$entercode==4,]))[1]    # 20 foreign installation
dim(unique(entrydat[entrydat$entercode==10,]))[1]   # 95 leader-reshuffling coup entry
dim(unique(entrydat[entrydat$entercode==11,]))[1]   # 88  (irregular) forced election/selection entry
dim(unique(entrydat[entrydat$entercode==12,]))[1]   # 172 first election/selection ("unprecedented" but regular/peaceful transition)
dim(unique(entrydat[entrydat$entercode==20,]))[1]  # 259 Regular selectino/election




##
## Exit manner

table(dat$exit, exclude=NULL )

# The variable "exit" codes exit manner for each "leader" (not leader-year). That is, the same exit manner is assigned to each leader and does not change during his entire tenure.
# If a leader still is in power in 2016, it is coded as "Still in Power".


#
#                  Assassination               Foreign Overthrow                     Leader Coup                   Legal Removal
#                            122                             202                             407                              59
#                  Natural Death                         Protest                     Regime Coup Resignation, Election/Selection
#                            841                             270                             501                            2138
#            Resignation, Health                  Still in Power            Successful Rebellion                            <NA>
#                             80                             968                             226                               0
#


exitdat <- dat[, c("ccode", "stateabb", "leadercode",  "leader", "exit", "eyear")    ]
# unique(exitdat)
dim(unique(exitdat))[1]    # 773 leaders (including puppet)

dim(unique(exitdat[exitdat$exit=="Assassination",]))[1]  # 19 assassination exit leaders

unique(exitdat[exitdat$exit=="Assassination",])
dim(unique(exitdat[exitdat$exit=="Foreign Overthrow",]))[1]   #24 foreign overthrow leaders
dim(unique(exitdat[exitdat$exit=="Leader Coup",]))[1]    # 83 reshuffling-coup exit leaders
dim(unique(exitdat[exitdat$exit=="Legal Removal",]))[1]   # 19 Legal removal leaders
dim(unique(exitdat[exitdat$exit=="Natural Death",]))[1]   # 53  natural death exit leaders
dim(unique(exitdat[exitdat$exit=="Protest",]))[1]   #  31 protest exit leaders
dim(unique(exitdat[exitdat$exit=="Regime Coup",]))[1]    # 66 regime-replacement coup exit leaders
dim(unique(exitdat[exitdat$exit=="Resignation, Election/Selection",]))[1]    # 378 voluntary resignation
dim(unique(exitdat[exitdat$exit=="Resignation, Health",]))[1]    # 8 leaders who resign due to health issue
dim(unique(exitdat[exitdat$exit=="Still in Power",]))[1]    # 64 leaders still in power (at the year 2016)
dim(unique(exitdat[exitdat$exit=="Successful Rebellion",]))[1]   # 28 successful rebel exit leaders



##
## Placeholder


# checking puppet

puppet <- dat[dat$puppet==1,]
length(table(puppet$leadercode))   # 93 puppet leaders

puppet <-  puppet[, c("ccode", "stateabb", "leadercode", "leader", "entry", "exit"  , "syear" , "eyear" ,  "puppet" )]
#unique(puppet) # list of all the puppet leaders
dim(unique(puppet))[1] # 93 puppets
















##
## Figure in the Manuscript
##

# We use only non-puppet leaders in our subsequent analyses

dat <- dat[dat$puppet!=1,]
dim(dat)    # [1] 5672   85


dat <- as.data.frame(dat)

attach(dat)


##
## Figure 1

setwd("C:/Users/pjb13215/Dropbox/Entry & Exit manner in Dictatorships")

betadraw <- read.table("Data/VCVmatrix_v1.csv", header=TRUE, sep=",")  # this is a 10000 draws of simulated coefficients.
dim(betadraw)

mnl1 <- multinom(multiexit  ~  reg_entry*log(tenure) + uncoord_entry*log(tenure) + foreign_entry*log(tenure)  +
 lngdp +growth + military + monarch + party + log_milper + age  + dwin + dlose + ddraw + incidencev414 , data=dat, Hess=TRUE)

mnl1.sum <- summary(mnl1)
mnl1.sum

coeffs1 <- mnl1.sum$coefficients
coeffs <- cbind(t(coeffs1[1,]), t(coeffs1[2,]), t(coeffs1[3,]), t(coeffs1[4,]), t(coeffs1[5,]))


# Set the tenure range for 1:10 years

tenure.range <- 1:10


# Irregular "Replacement" Entry leader
Xhyp   <- cbind(1,   0, log(tenure.range), 0, 0,  mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE)  ,median(monarch, na.rm=TRUE) ,median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 0*log(tenure.range),  0*log(tenure.range),   0*log(tenure.range))

dim(Xhyp)


bgr2 <- Xhyp%*%coeffs[1:19]
bgr3 <- Xhyp%*%coeffs[20:38]
bgr4 <- Xhyp%*%coeffs[39:57]
bgr5 <- Xhyp%*%coeffs[58:76]
bgr6 <- Xhyp%*%coeffs[77:95]


prob1gr <- exp(0) / (exp(0) + exp(bgr2) + exp(bgr3) + exp(bgr4) + exp(bgr5) + exp(bgr6))
prob2gr <- exp(bgr2) /(exp(0) + exp(bgr2) + exp(bgr3) + exp(bgr4) + exp(bgr5) + exp(bgr6))
prob3gr <- exp(bgr3) / (exp(0) + exp(bgr2) + exp(bgr3) + exp(bgr4) + exp(bgr5) + exp(bgr6))
prob4gr <- exp(bgr4) / (exp(0) + exp(bgr2) + exp(bgr3) + exp(bgr4) + exp(bgr5) + exp(bgr6))
prob5gr <- exp(bgr5) /(exp(0) + exp(bgr2) + exp(bgr3) + exp(bgr4) + exp(bgr5) + exp(bgr6))
prob6gr <- exp(bgr6) /  (exp(0) + exp(bgr2) + exp(bgr3) + exp(bgr4) + exp(bgr5) + exp(bgr6))



xbgr2.sim <- Xhyp%*% t(betadraw[,1:19])
xbgr3.sim <- Xhyp%*% t(betadraw[,20:38])
xbgr4.sim <- Xhyp%*% t(betadraw[,39:57])
xbgr5.sim <- Xhyp%*% t(betadraw[,58:76])
xbgr6.sim <- Xhyp%*% t(betadraw[,77:95])


prob2gr.sim <- exp(xbgr2.sim)/ (exp(0) + exp(xbgr2.sim) + exp(xbgr3.sim) + exp(xbgr4.sim) + exp(xbgr5.sim) + exp(xbgr6.sim) )
prob3gr.sim <- exp(xbgr3.sim)/ (exp(0) + exp(xbgr2.sim) + exp(xbgr3.sim) + exp(xbgr4.sim) + exp(xbgr5.sim) + exp(xbgr6.sim) )
prob4gr.sim <- exp(xbgr4.sim)/ (exp(0) + exp(xbgr2.sim) + exp(xbgr3.sim) + exp(xbgr4.sim) + exp(xbgr5.sim) + exp(xbgr6.sim) )
prob5gr.sim <- exp(xbgr5.sim)/ (exp(0) + exp(xbgr2.sim) + exp(xbgr3.sim) + exp(xbgr4.sim) + exp(xbgr5.sim) + exp(xbgr6.sim) )
prob6gr.sim <- exp(xbgr6.sim)/ (exp(0) + exp(xbgr2.sim) + exp(xbgr3.sim) + exp(xbgr4.sim) + exp(xbgr5.sim) + exp(xbgr6.sim) )



# Irregular "Reorganization" Entry Leader
Yhyp   <- cbind(1,   0, log(tenure.range), 1, 0,  mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE)  ,median(monarch, na.rm=TRUE) ,median(party, na.rm=TRUE) , mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 0*log(tenure.range),  1*log(tenure.range),   0*log(tenure.range))


ybgr2 <- Yhyp%*%coeffs[1:19]
ybgr3 <- Yhyp%*%coeffs[20:38]
ybgr4 <- Yhyp%*%coeffs[39:57]
ybgr5 <- Yhyp%*%coeffs[58:76]
ybgr6 <- Yhyp%*%coeffs[77:95]

prob1ygr <- exp(0) / (exp(0) + exp(ybgr2) + exp(ybgr3) + exp(ybgr4)+ exp(ybgr5) + exp(ybgr6))
prob2ygr <- exp(ybgr2) /(exp(0) + exp(ybgr2) + exp(ybgr3) + exp(ybgr4)+ exp(ybgr5) + exp(ybgr6))
prob3ygr <- exp(ybgr3) /(exp(0) + exp(ybgr2) + exp(ybgr3) + exp(ybgr4)+ exp(ybgr5) + exp(ybgr6))
prob4ygr <- exp(ybgr4) /(exp(0) + exp(ybgr2) + exp(ybgr3) + exp(ybgr4)+ exp(ybgr5) + exp(ybgr6))
prob5ygr <- exp(ybgr5) /(exp(0) + exp(ybgr2) + exp(ybgr3) + exp(ybgr4)+ exp(ybgr5) + exp(ybgr6))
prob6ygr <- exp(ybgr6) /(exp(0) + exp(ybgr2) + exp(ybgr3) + exp(ybgr4)+ exp(ybgr5) + exp(ybgr6))


ybgr2.sim <- Yhyp%*% t(betadraw[,1:19])
ybgr3.sim <- Yhyp%*% t(betadraw[,20:38])
ybgr4.sim <- Yhyp%*% t(betadraw[,39:57])
ybgr5.sim <- Yhyp%*% t(betadraw[,58:76])
ybgr6.sim <- Yhyp%*% t(betadraw[,77:95])


prob2ygr.sim <- exp(ybgr2.sim )/ (exp(0) + exp(ybgr2.sim) + exp(ybgr3.sim ) + exp(ybgr4.sim)+ exp(ybgr5.sim) + exp(ybgr6.sim)  )
prob3ygr.sim <- exp(ybgr3.sim )/ (exp(0) + exp(ybgr2.sim) + exp(ybgr3.sim ) + exp(ybgr4.sim)+ exp(ybgr5.sim) + exp(ybgr6.sim)  )
prob4ygr.sim <- exp(ybgr4.sim )/ (exp(0) + exp(ybgr2.sim) + exp(ybgr3.sim ) + exp(ybgr4.sim)+ exp(ybgr5.sim) + exp(ybgr6.sim)  )
prob5ygr.sim <- exp(ybgr5.sim )/ (exp(0) + exp(ybgr2.sim) + exp(ybgr3.sim ) + exp(ybgr4.sim)+ exp(ybgr5.sim) + exp(ybgr6.sim)  )
prob6ygr.sim <- exp(ybgr6.sim )/ (exp(0) + exp(ybgr2.sim) + exp(ybgr3.sim ) + exp(ybgr4.sim)+ exp(ybgr5.sim) + exp(ybgr6.sim)  )




# Regular Entry leader
Zhyp <-  cbind(1,   1, log(tenure.range), 0, 0,  mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE)  ,median(monarch, na.rm=TRUE) ,median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 1*log(tenure.range),  0*log(tenure.range),   0*log(tenure.range))


zbgr2 <- Zhyp%*%coeffs[1:19]
zbgr3 <- Zhyp%*%coeffs[20:38]
zbgr4 <- Zhyp%*%coeffs[39:57]
zbgr5 <- Zhyp%*%coeffs[58:76]
zbgr6 <- Zhyp%*%coeffs[77:95]


prob1zgr <- exp(0) / (exp(0) + exp(zbgr2) + exp(zbgr3) + exp(zbgr4)+ exp(zbgr5) + exp(zbgr6))
prob2zgr <- exp(zbgr2) / (exp(0) + exp(zbgr2) + exp(zbgr3) + exp(zbgr4)+ exp(zbgr5) + exp(zbgr6))
prob3zgr <- exp(zbgr3) / (exp(0) + exp(zbgr2) + exp(zbgr3) + exp(zbgr4)+ exp(zbgr5) + exp(zbgr6))
prob4zgr <- exp(zbgr4) / (exp(0) + exp(zbgr2) + exp(zbgr3) + exp(zbgr4)+ exp(zbgr5) + exp(zbgr6))
prob5zgr <- exp(zbgr5) / (exp(0) + exp(zbgr2) + exp(zbgr3) + exp(zbgr4)+ exp(zbgr5) + exp(zbgr6))
prob6zgr <- exp(zbgr6) / (exp(0) + exp(zbgr2) + exp(zbgr3) + exp(zbgr4)+ exp(zbgr5) + exp(zbgr6))



zbgr2.sim <- Zhyp%*% t(betadraw[,1:19])
zbgr3.sim <- Zhyp%*% t(betadraw[,20:38])
zbgr4.sim <- Zhyp%*% t(betadraw[,39:57])
zbgr5.sim <- Zhyp%*% t(betadraw[,58:76])
zbgr6.sim <- Zhyp%*% t(betadraw[,77:95])

prob2zgr.sim <- exp(zbgr2.sim )/ (exp(0) + exp(zbgr2.sim) + exp(zbgr3.sim ) + exp(zbgr4.sim) + exp(zbgr5.sim ) + exp(zbgr6.sim)  )
prob3zgr.sim <- exp(zbgr3.sim )/ (exp(0) + exp(zbgr2.sim) + exp(zbgr3.sim ) + exp(zbgr4.sim) + exp(zbgr5.sim ) + exp(zbgr6.sim)  )
prob4zgr.sim <- exp(zbgr4.sim )/ (exp(0) + exp(zbgr2.sim) + exp(zbgr3.sim ) + exp(zbgr4.sim) + exp(zbgr5.sim ) + exp(zbgr6.sim)  )
prob5zgr.sim <- exp(zbgr5.sim )/ (exp(0) + exp(zbgr2.sim) + exp(zbgr3.sim ) + exp(zbgr4.sim) + exp(zbgr5.sim ) + exp(zbgr6.sim)  )
prob6zgr.sim <- exp(zbgr6.sim )/ (exp(0) + exp(zbgr2.sim) + exp(zbgr3.sim ) + exp(zbgr4.sim) + exp(zbgr5.sim ) + exp(zbgr6.sim)  )




## First Difference :

# Pr(Coalition-Competing Removal|Replacement entry) - Pr(Coalition-Competing Removal|Regular entry)
Diff.P.sim.gy <- prob2gr.sim - prob2ygr.sim

Diff.P.gy <- prob2gr - prob2ygr

# Pr(Coalition-Competing Removal|Replacement Entry ) - Pr(Coalition-Competing Exit| Irregular Reorganization Entry)
Diff.P.sim.gz <- prob2gr.sim -  prob2zgr.sim
Diff.P.gz <-  prob2gr  -  prob2zgr



##
## Figure 1
##

windows()
par(mar=c(4,4,3,1), mfrow=c(1,2))
# A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot.
# The default is c(5, 4, 4, 2) + 0.1.


plot(tenure.range, rowMeans(Diff.P.sim.gz) ,type="l", xlim=c(1, 10),  xaxt='n', ylim=c(-0.015 , 0.02),     lwd=3,ylab="Pr(Competing Exit|Replacement Entry) - Pr(Competing Exit|Regular Entry)", xlab="Tenure")
lines(tenure.range, t(apply(Diff.P.sim.gz,1,sort))[,9750],type="l",lty=3,lwd=2    )
lines(tenure.range,  t(apply(Diff.P.sim.gz,1,sort))[,250],type="l",lty=3,lwd=2  )
abline(h=0,lty=3 )
axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T,  line=-2)
legend(1, 0.02,  bty="n",   c( "First Difference (v.s. Regular Entry)" , "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")



plot(tenure.range, rowMeans(Diff.P.sim.gy) ,type="l", xlim=c(1, 10), xaxt='n', ylim=c(-0.03 , 0.02),   lwd=3,ylab="Pr(Competing Exit|Replacement Entry) - Pr(Competing Exit|Reorganization Entry)", xlab="Tenure")
lines(tenure.range, t(apply(Diff.P.sim.gy,1,sort))[,9750],type="l",lty=3,lwd=2   )
lines(tenure.range,  t(apply(Diff.P.sim.gy,1,sort))[,250],type="l",lty=3,lwd=2  )
#legend(x=500, y=0.6, bty="n", legend="Incumbency Advantage", lty=1, lwd=3, seg.len=0.75, col="blue", x.intersp=0.5)
abline(h=0,lty=3 )

axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)
legend(1, 0.02,  bty="n",   c( "First Difference (v.s. Reorganization Entry)" , "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")



# save

savePlot(filename="ISQ Submission/R&R/meffect_competing_removal", type = "pdf")









## First Difference of Other Irregular Exit (Coalition-collapse & Circumventing exit):
# Pr(Other Irregular exit|Replacement entry) - Pr(Other Irregular Exit|Reorganization entry)

Diff.P3.sim.gy <- prob3gr.sim - prob3ygr.sim

Diff.P4.sim.gy <- prob4gr.sim - prob4ygr.sim

# First-Difference b/w Coordinated Irregular Entry and Regular Entry


Diff.P3.sim.gz <-   prob3gr.sim - prob3zgr.sim


Diff.P4.sim.gz <-   prob4gr.sim - prob4zgr.sim





##
## Figure 2
##

# DV: Coalition-Circumventing Exit & Coalition-Collapsing Exit (Hypothesis 2 )


windows()
par(mar=c(4,4,3,1), mfrow=c(2,2))
# A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot.
# The default is c(5, 4, 4, 2) + 0.1.


plot(tenure.range, rowMeans(Diff.P3.sim.gz) ,type="l", xlim=c(1, 10),main="DV: Circumventing Exit", xaxt='n', ylim=c(-0.01, 0.04), lwd=3,ylab="First Difference", xlab="Tenure")
lines(tenure.range, t(apply(Diff.P3.sim.gz,1,sort))[,9750],type="l",lty=3,lwd=2    )
lines(tenure.range,  t(apply(Diff.P3.sim.gz,1,sort))[,250],type="l",lty=3,lwd=2  )
abline(h=0,lty=3 )
axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T,  line=-2)
legend(1, 0.04,  bty="n",   c( "First Difference (v.s. Regular Entry)" , "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")



plot(tenure.range, rowMeans(Diff.P3.sim.gy) ,type="l", xlim=c(1, 10),main="DV: Circumventing Exit", xaxt='n', ylim=c(-0.06 , 0.04),   lwd=3,ylab="First Difference", xlab="Tenure")
lines(tenure.range, t(apply(Diff.P3.sim.gy,1,sort))[,9750],type="l",lty=3,lwd=2   )
lines(tenure.range,  t(apply(Diff.P3.sim.gy,1,sort))[,250],type="l",lty=3,lwd=2  )
#legend(x=500, y=0.6, bty="n", legend="Incumbency Advantage", lty=1, lwd=3, seg.len=0.75, col="blue", x.intersp=0.5)
abline(h=0,lty=3 )

axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)
legend(1, 0.04,  bty="n",   c( "First Difference (v.s. Reorganization Entry)" , "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")




plot(tenure.range, rowMeans(Diff.P4.sim.gz) ,type="l", xlim=c(1, 10),main="DV: Collapsing Exit", xaxt='n', ylim=c(-0.03, 0.03), lwd=3,ylab="First Difference", xlab="Tenure")
lines(tenure.range, t(apply(Diff.P4.sim.gz,1,sort))[,9750],type="l",lty=3,lwd=2    )
lines(tenure.range,  t(apply(Diff.P4.sim.gz,1,sort))[,250],type="l",lty=3,lwd=2  )
abline(h=0,lty=3 )
axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T,  line=-2)
legend(1, 0.03,  bty="n",   c( "First Difference (v.s. Regular Entry)" , "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")



plot(tenure.range, rowMeans(Diff.P4.sim.gy) ,type="l", xlim=c(1, 10),main="DV: Collapsing Exit", xaxt='n', ylim=c(-0.05 , 0.03),   lwd=3,ylab="First Difference", xlab="Tenure")
lines(tenure.range, t(apply(Diff.P4.sim.gy,1,sort))[,9750],type="l",lty=3,lwd=2   )
lines(tenure.range,  t(apply(Diff.P4.sim.gy,1,sort))[,250],type="l",lty=3,lwd=2  )
#legend(x=500, y=0.6, bty="n", legend="Incumbency Advantage", lty=1, lwd=3, seg.len=0.75, col="blue", x.intersp=0.5)
abline(h=0,lty=3 )

axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)
legend(1, 0.03,  bty="n",   c( "First Difference (v.s. Reorganization Entry)" , "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")





# save

savePlot(filename="ISQ Submission/R&R/meffect_circumventing_removal", type = "pdf")















##
## Figure 3
##

# Collapsing IV ( Entry Manner):

# Model 2 (Base; regular entry) :

betadraw <- read.table("C:/Users/pjb13215/Dropbox/Entry & Exit manner in Dictatorships/Data/VCVmatrix_v2.csv", header=TRUE, sep=",")  # this is a 10000 draws of simulated coefficients.
dim(betadraw) # 10000 * 12

mnl1 <- multinom(multiexit ~  rebelentry*log(tenure) + regimecoup_entry*log(tenure) + leadercoup_entry*log(tenure) + ireg_election_entry*log(tenure) +
foreign_entry*log(tenure)  + lngdp +growth + military + monarch + party + log_milper + age  + dwin + dlose + ddraw + incidencev414 , data=dat, Hess=TRUE)

mnl1.sum <- summary(mnl1)
mnl1.sum

coeffs1 <- mnl1.sum$coefficients
coeffs <- cbind(t(coeffs1[1,]), t(coeffs1[2,]), t(coeffs1[3,]), t(coeffs1[4,]), t(coeffs1[5,]))

# Set the tenure range (Figure)

tenure.range <- 1:10

# Rebel Entry leader
Xhyp   <- cbind(1,   1, log(tenure.range), 0, 0, 0 ,0,   mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 1*log(tenure.range),  0*log(tenure.range),   0*log(tenure.range), 0*log(tenure.range),   0*log(tenure.range))


dim(Xhyp)

bgr2 <- Xhyp%*%coeffs[1:23]
bgr3 <- Xhyp%*%coeffs[24:46]
bgr4 <- Xhyp%*%coeffs[47:69]
bgr5 <- Xhyp%*%coeffs[70:92]
bgr6 <- Xhyp%*%coeffs[93:115]

prob1gr <- exp(0) / (exp(0) + exp(bgr2) + exp(bgr3) + exp(bgr4) + exp(bgr5)+ exp(bgr6)  )
prob2gr <- exp(bgr2) /  (exp(0) + exp(bgr2) + exp(bgr3) + exp(bgr4) + exp(bgr5)+ exp(bgr6)  )
prob3gr <- exp(bgr3) /  (exp(0) + exp(bgr2) + exp(bgr3) + exp(bgr4) + exp(bgr5)+ exp(bgr6)  )
prob4gr <- exp(bgr4) /   (exp(0) + exp(bgr2) + exp(bgr3) + exp(bgr4) + exp(bgr5)+ exp(bgr6)  )
prob5gr <- exp(bgr5) /   (exp(0) + exp(bgr2) + exp(bgr3) + exp(bgr4) + exp(bgr5)+ exp(bgr6)  )
prob6gr <- exp(bgr6) /   (exp(0) + exp(bgr2) + exp(bgr3) + exp(bgr4) + exp(bgr5)+ exp(bgr6)  )


xbgr2.sim <- Xhyp%*% t(betadraw[,1:23])
xbgr3.sim <- Xhyp%*% t(betadraw[,24:46])
xbgr4.sim <- Xhyp%*% t(betadraw[,47:69])
xbgr5.sim <- Xhyp%*% t(betadraw[,70:92])
xbgr6.sim <- Xhyp%*% t(betadraw[,93:115])


prob2gr.sim <- exp(xbgr2.sim)/ (exp(0) + exp(xbgr2.sim) + exp(xbgr3.sim) + exp(xbgr4.sim) + exp(xbgr5.sim) + exp(xbgr6.sim) )
prob3gr.sim <- exp(xbgr3.sim)/ (exp(0) + exp(xbgr2.sim) + exp(xbgr3.sim) + exp(xbgr4.sim) + exp(xbgr5.sim) + exp(xbgr6.sim) )
prob4gr.sim <- exp(xbgr4.sim)/ (exp(0) + exp(xbgr2.sim) + exp(xbgr3.sim) + exp(xbgr4.sim) + exp(xbgr5.sim) + exp(xbgr6.sim) )
prob5gr.sim <- exp(xbgr5.sim)/ (exp(0) + exp(xbgr2.sim) + exp(xbgr3.sim) + exp(xbgr4.sim) + exp(xbgr5.sim) + exp(xbgr6.sim) )
prob6gr.sim <- exp(xbgr6.sim)/ (exp(0) + exp(xbgr2.sim) + exp(xbgr3.sim) + exp(xbgr4.sim) + exp(xbgr5.sim) + exp(xbgr6.sim) )




# Regime-Coup Entry Leader
Yhyp  <- cbind(1,   0, log(tenure.range), 1, 0, 0 ,0,   mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 0*log(tenure.range),  1*log(tenure.range),   0*log(tenure.range), 0*log(tenure.range),   0*log(tenure.range))

ybgr2 <- Yhyp%*%coeffs[1:23]
ybgr3 <- Yhyp%*%coeffs[24:46]
ybgr4 <- Yhyp%*%coeffs[47:69]
ybgr5 <- Yhyp%*%coeffs[70:92]
ybgr6 <- Yhyp%*%coeffs[93:115]


prob1ygr <- exp(0) / (exp(0) + exp(ybgr2) + exp(ybgr3) + exp(ybgr4)+ exp(ybgr5) + exp(ybgr6))
prob2ygr <- exp(ybgr2) /(exp(0) + exp(ybgr2) + exp(ybgr3) + exp(ybgr4)+ exp(ybgr5) + exp(ybgr6))
prob3ygr <- exp(ybgr3) /(exp(0) + exp(ybgr2) + exp(ybgr3) + exp(ybgr4)+ exp(ybgr5) + exp(ybgr6))
prob4ygr <- exp(ybgr4) /(exp(0) + exp(ybgr2) + exp(ybgr3) + exp(ybgr4)+ exp(ybgr5) + exp(ybgr6))
prob5ygr <- exp(ybgr5) /(exp(0) + exp(ybgr2) + exp(ybgr3) + exp(ybgr4)+ exp(ybgr5) + exp(ybgr6))
prob6ygr <- exp(ybgr6) /(exp(0) + exp(ybgr2) + exp(ybgr3) + exp(ybgr4)+ exp(ybgr5) + exp(ybgr6))



ybgr2.sim <- Yhyp%*% t(betadraw[,1:23])
ybgr3.sim <- Yhyp%*% t(betadraw[,24:46])
ybgr4.sim <- Yhyp%*% t(betadraw[,47:69])
ybgr5.sim <- Yhyp%*% t(betadraw[,70:92])
ybgr6.sim <- Yhyp%*% t(betadraw[,93:115])


prob2ygr.sim <- exp(ybgr2.sim)/ (exp(0) + exp(ybgr2.sim) + exp(ybgr3.sim) + exp(ybgr4.sim) + exp(ybgr5.sim) + exp(ybgr6.sim) )
prob3ygr.sim <- exp(ybgr3.sim)/ (exp(0) + exp(ybgr2.sim) + exp(ybgr3.sim) + exp(ybgr4.sim) + exp(ybgr5.sim) + exp(ybgr6.sim) )
prob4ygr.sim <- exp(ybgr4.sim)/ (exp(0) + exp(ybgr2.sim) + exp(ybgr3.sim) + exp(ybgr4.sim) + exp(ybgr5.sim) + exp(ybgr6.sim) )
prob5ygr.sim <- exp(ybgr5.sim)/ (exp(0) + exp(ybgr2.sim) + exp(ybgr3.sim) + exp(ybgr4.sim) + exp(ybgr5.sim) + exp(ybgr6.sim) )
prob6ygr.sim <- exp(ybgr6.sim)/(exp(0) + exp(ybgr2.sim) + exp(ybgr3.sim) + exp(ybgr4.sim) + exp(ybgr5.sim) + exp(ybgr6.sim) )


# Regular leader
Zhyp <- cbind(1,   0, log(tenure.range), 0, 0, 0 ,0,   mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 0*log(tenure.range),  0*log(tenure.range),   0*log(tenure.range), 0*log(tenure.range),   0*log(tenure.range))


zbgr2 <- Zhyp%*%coeffs[1:23]
zbgr3 <- Zhyp%*%coeffs[24:46]
zbgr4 <- Zhyp%*%coeffs[47:69]
zbgr5 <- Zhyp%*%coeffs[70:92]
zbgr6 <- Zhyp%*%coeffs[93:115]


prob1zgr <- exp(0) / (exp(0) + exp(zbgr2) + exp(zbgr3) + exp(zbgr4)+ exp(zbgr5) + exp(zbgr6))
prob2zgr <- exp(zbgr2) / (exp(0) + exp(zbgr2) + exp(zbgr3) + exp(zbgr4)+ exp(zbgr5) + exp(zbgr6))
prob3zgr <- exp(zbgr3) / (exp(0) + exp(zbgr2) + exp(zbgr3) + exp(zbgr4)+ exp(zbgr5) + exp(zbgr6))
prob4zgr <- exp(zbgr4) /(exp(0) + exp(zbgr2) + exp(zbgr3) + exp(zbgr4)+ exp(zbgr5) + exp(zbgr6))
prob5zgr <- exp(zbgr5) /(exp(0) + exp(zbgr2) + exp(zbgr3) + exp(zbgr4)+ exp(zbgr5) + exp(zbgr6))
prob6zgr <- exp(zbgr6) /(exp(0) + exp(zbgr2) + exp(zbgr3) + exp(zbgr4)+ exp(zbgr5) + exp(zbgr6))


zbgr2.sim <- Zhyp%*% t(betadraw[,1:23])
zbgr3.sim <- Zhyp%*% t(betadraw[,24:46])
zbgr4.sim <- Zhyp%*% t(betadraw[,47:69])
zbgr5.sim <- Zhyp%*% t(betadraw[,70:92])
zbgr6.sim <- Zhyp%*% t(betadraw[,93:115])


prob2zgr.sim <- exp(zbgr2.sim)/ (exp(0) + exp(zbgr2.sim) + exp(zbgr3.sim) + exp(zbgr4.sim) + exp(zbgr5.sim) + exp(zbgr6.sim) )
prob3zgr.sim <- exp(zbgr3.sim)/  (exp(0) + exp(zbgr2.sim) + exp(zbgr3.sim) + exp(zbgr4.sim) + exp(zbgr5.sim) + exp(zbgr6.sim) )
prob4zgr.sim <- exp(zbgr4.sim)/  (exp(0) + exp(zbgr2.sim) + exp(zbgr3.sim) + exp(zbgr4.sim) + exp(zbgr5.sim) + exp(zbgr6.sim) )
prob5zgr.sim <- exp(zbgr5.sim)/  (exp(0) + exp(zbgr2.sim) + exp(zbgr3.sim) + exp(zbgr4.sim) + exp(zbgr5.sim) + exp(zbgr6.sim) )
prob6zgr.sim <- exp(zbgr6.sim)/ (exp(0) + exp(zbgr2.sim) + exp(zbgr3.sim) + exp(zbgr4.sim) + exp(zbgr5.sim) + exp(zbgr6.sim) )



# Leader-Reshuffle Coup Leader
Rhyp  <- cbind(1,   0, log(tenure.range), 0, 1, 0 ,0,   mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 0*log(tenure.range),  0*log(tenure.range),   1*log(tenure.range), 0*log(tenure.range),   0*log(tenure.range))


rbgr2 <- Rhyp%*%coeffs[1:23]
rbgr3 <- Rhyp%*%coeffs[24:46]
rbgr4 <- Rhyp%*%coeffs[47:69]
rbgr5 <- Rhyp%*%coeffs[70:92]
rbgr6 <- Rhyp%*%coeffs[93:115]


prob1rgr <- exp(0) / (exp(0) + exp(rbgr2) + exp(rbgr3) + exp(rbgr4)+ exp(rbgr5) + exp(rbgr6) )
prob2rgr <- exp(rbgr2) /  (exp(0) + exp(rbgr2) + exp(rbgr3) + exp(rbgr4)+ exp(rbgr5) + exp(rbgr6) )
prob3rgr <- exp(rbgr3) /  (exp(0) + exp(rbgr2) + exp(rbgr3) + exp(rbgr4)+ exp(rbgr5) + exp(rbgr6) )
prob4rgr <- exp(rbgr4) / (exp(0) + exp(rbgr2) + exp(rbgr3) + exp(rbgr4)+ exp(rbgr5) + exp(rbgr6) )
prob5rgr <- exp(rbgr5) / (exp(0) + exp(rbgr2) + exp(rbgr3) + exp(rbgr4)+ exp(rbgr5) + exp(rbgr6) )
prob6rgr <- exp(rbgr6) / (exp(0) + exp(rbgr2) + exp(rbgr3) + exp(rbgr4)+ exp(rbgr5) + exp(rbgr6) )


rbgr2.sim <- Rhyp%*% t(betadraw[,1:23])
rbgr3.sim <- Rhyp%*% t(betadraw[,24:46])
rbgr4.sim <- Rhyp%*% t(betadraw[,47:69])
rbgr5.sim <- Rhyp%*% t(betadraw[,70:92])
rbgr6.sim <- Rhyp%*% t(betadraw[,93:115])


prob2rgr.sim <- exp(rbgr2.sim)/ (exp(0) + exp(rbgr2.sim) + exp(rbgr3.sim) + exp(rbgr4.sim) + exp(rbgr5.sim) + exp(rbgr6.sim) )
prob3rgr.sim <- exp(rbgr3.sim)/ (exp(0) + exp(rbgr2.sim) + exp(rbgr3.sim) + exp(rbgr4.sim) + exp(rbgr5.sim) + exp(rbgr6.sim) )
prob4rgr.sim <- exp(rbgr4.sim)/ (exp(0) + exp(rbgr2.sim) + exp(rbgr3.sim) + exp(rbgr4.sim) + exp(rbgr5.sim) + exp(rbgr6.sim) )
prob5rgr.sim <- exp(rbgr5.sim)/  (exp(0) + exp(rbgr2.sim) + exp(rbgr3.sim) + exp(rbgr4.sim) + exp(rbgr5.sim) + exp(rbgr6.sim) )
prob6rgr.sim <- exp(rbgr6.sim)/  (exp(0) + exp(rbgr2.sim) + exp(rbgr3.sim) + exp(rbgr4.sim) + exp(rbgr5.sim) + exp(rbgr6.sim) )






# Other Irregular Leader
Ohyp <- cbind(1,   0, log(tenure.range), 0, 0, 1 ,0,   mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 0*log(tenure.range),  0*log(tenure.range),   0*log(tenure.range), 1*log(tenure.range),   0*log(tenure.range))


obgr2 <- Ohyp%*%coeffs[1:23]
obgr3 <- Ohyp%*%coeffs[24:46]
obgr4 <- Ohyp%*%coeffs[47:69]
obgr5 <- Ohyp%*%coeffs[70:92]
obgr6 <- Ohyp%*%coeffs[93:115]

prob1ogr <- exp(0) / (exp(0) + exp(obgr2) + exp(obgr3) + exp(obgr4)+ exp(obgr5) + exp(obgr6))
prob2ogr <- exp(obgr2) / (exp(0) + exp(obgr2) + exp(obgr3) + exp(obgr4)+ exp(obgr5) + exp(obgr6))
prob3ogr <- exp(obgr3) /(exp(0) + exp(obgr2) + exp(obgr3) + exp(obgr4)+ exp(obgr5) + exp(obgr6))
prob4ogr <- exp(obgr4) /(exp(0) + exp(obgr2) + exp(obgr3) + exp(obgr4)+ exp(obgr5) + exp(obgr6))
prob5ogr <- exp(obgr5) /(exp(0) + exp(obgr2) + exp(obgr3) + exp(obgr4)+ exp(obgr5) + exp(obgr6))
prob6ogr <- exp(obgr6) /(exp(0) + exp(obgr2) + exp(obgr3) + exp(obgr4)+ exp(obgr5) + exp(obgr6))


obgr2.sim <- Ohyp%*% t(betadraw[,1:23])
obgr3.sim <- Ohyp%*% t(betadraw[,24:46])
obgr4.sim <- Ohyp%*% t(betadraw[,47:69])
obgr5.sim <- Ohyp%*% t(betadraw[,70:92])
obgr6.sim <- Ohyp%*% t(betadraw[,93:115])


prob2ogr.sim <- exp(obgr2.sim)/ (exp(0) + exp(obgr2.sim) + exp(obgr3.sim) + exp(obgr4.sim) + exp(obgr5.sim) + exp(obgr6.sim) )
prob3ogr.sim <- exp(obgr3.sim)/ (exp(0) + exp(obgr2.sim) + exp(obgr3.sim) + exp(obgr4.sim) + exp(obgr5.sim) + exp(obgr6.sim) )
prob4ogr.sim <- exp(obgr4.sim)/(exp(0) + exp(obgr2.sim) + exp(obgr3.sim) + exp(obgr4.sim) + exp(obgr5.sim) + exp(obgr6.sim) )
prob5ogr.sim <- exp(obgr5.sim)/ (exp(0) + exp(obgr2.sim) + exp(obgr3.sim) + exp(obgr4.sim) + exp(obgr5.sim) + exp(obgr6.sim) )
prob6ogr.sim <- exp(obgr6.sim)/ (exp(0) + exp(obgr2.sim) + exp(obgr3.sim) + exp(obgr4.sim) + exp(obgr5.sim) + exp(obgr6.sim) )




## Coalition-Competing Exit
# Pr(Coalition Competing|Rebel entry) - Pr(Coalition Competing|Regular entry)

Diff.P2.sim.gz <-  prob2gr.sim - prob2zgr.sim

# Pr(Coalition Competing|Regime-Change entry) - Pr(Coalition Competing|Regular entry)

Diff.P2.sim.yz <-  prob2ygr.sim - prob2zgr.sim


# Pr(Coalition Competing|Reshuffling coup entry) - Pr(Coalition Competing|Regular entry)

Diff.P2.sim.rz <- prob2rgr.sim - prob2zgr.sim

# Pr(Coalition Competing|Force Election entry) - Pr(Coalition Competing|Regular entry)

Diff.P2.sim.oz <- prob2ogr.sim - prob2zgr.sim




##
## Model 3 (base category :Irregular Reorganization Entry) :
##

betadraw3 <- read.table("C:/Users/pjb13215/Dropbox/Entry & Exit manner in Dictatorships/Data/VCVmatrix_v2_2.csv", header=TRUE, sep=",")  # this is a 10000 draws of simulated coefficients.
dim(betadraw3) # 10000 * 12

mnl3 <- multinom(multiexit ~  rebelentry*log(tenure) + regimecoup_entry*log(tenure) +
foreign_entry*log(tenure)  + first_election_entry*log(tenure) + reg_election_entry *log(tenure) +
lngdp +growth + military + monarch + party + log_milper + age  + dwin + dlose + ddraw + incidencev414 , data=dat, Hess=TRUE)


mnl3.sum <- summary(mnl3)
mnl3.sum

coeffs3 <- mnl3.sum$coefficients
coeffs3 <- cbind(t(coeffs3[1,]), t(coeffs3[2,]), t(coeffs3[3,]), t(coeffs3[4,]), t(coeffs3[5,]))
#covmat <- solve(mnl1$Hessian)

# Set the tenure range (Figure)

tenure.range <- 1:10

# Rebel Entry leader
Ahyp   <- cbind(1,   1, log(tenure.range), 0, 0, 0 ,0,   mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 1*log(tenure.range),  0*log(tenure.range),   0*log(tenure.range), 0*log(tenure.range),   0*log(tenure.range))


abgr2.sim <- Ahyp%*% t(betadraw3[,1:23])
abgr3.sim <- Ahyp%*% t(betadraw3[,24:46])
abgr4.sim <- Ahyp%*% t(betadraw3[,47:69])
abgr5.sim <- Ahyp%*% t(betadraw3[,70:92])
abgr6.sim <- Ahyp%*% t(betadraw3[,93:115])


prob2agr.sim <- exp(abgr2.sim)/ (exp(0) + exp(abgr2.sim) + exp(abgr3.sim) + exp(abgr4.sim) + exp(abgr5.sim) + exp(abgr6.sim) )
prob3agr.sim <- exp(abgr3.sim)/(exp(0) + exp(abgr2.sim) + exp(abgr3.sim) + exp(abgr4.sim) + exp(abgr5.sim) + exp(abgr6.sim) )
prob4agr.sim <- exp(abgr4.sim)/ (exp(0) + exp(abgr2.sim) + exp(abgr3.sim) + exp(abgr4.sim) + exp(abgr5.sim) + exp(abgr6.sim) )
prob5agr.sim <- exp(abgr5.sim)/(exp(0) + exp(abgr2.sim) + exp(abgr3.sim) + exp(abgr4.sim) + exp(abgr5.sim) + exp(abgr6.sim) )
prob6agr.sim <- exp(abgr6.sim)/(exp(0) + exp(abgr2.sim) + exp(abgr3.sim) + exp(abgr4.sim) + exp(abgr5.sim) + exp(abgr6.sim) )



# Regime-Coup Entry Leader
Bhyp   <- cbind(1,   0, log(tenure.range), 1, 0, 0 ,0,   mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) , 0*log(tenure.range),  1*log(tenure.range),   0*log(tenure.range), 0*log(tenure.range),   0*log(tenure.range))


bbgr2.sim <- Bhyp%*% t(betadraw3[,1:23])
bbgr3.sim <- Bhyp%*% t(betadraw3[,24:46])
bbgr4.sim <- Bhyp%*% t(betadraw3[,47:69])
bbgr5.sim <- Bhyp%*% t(betadraw3[,70:92])
bbgr6.sim <- Bhyp%*% t(betadraw3[,93:115])


prob2bgr.sim <- exp(bbgr2.sim)/ (exp(0) + exp(bbgr2.sim) + exp(bbgr3.sim) + exp(bbgr4.sim) + exp(bbgr5.sim) + exp(bbgr6.sim) )
prob3bgr.sim <- exp(bbgr3.sim)/(exp(0) + exp(bbgr2.sim) + exp(bbgr3.sim) + exp(bbgr4.sim) + exp(bbgr5.sim) + exp(bbgr6.sim) )
prob4bgr.sim <- exp(bbgr4.sim)/(exp(0) + exp(bbgr2.sim) + exp(bbgr3.sim) + exp(bbgr4.sim) + exp(bbgr5.sim) + exp(bbgr6.sim) )
prob5bgr.sim <- exp(bbgr5.sim)/(exp(0) + exp(bbgr2.sim) + exp(bbgr3.sim) + exp(bbgr4.sim) + exp(bbgr5.sim) + exp(bbgr6.sim) )
prob6bgr.sim <- exp(bbgr6.sim)/(exp(0) + exp(bbgr2.sim) + exp(bbgr3.sim) + exp(bbgr4.sim) + exp(bbgr5.sim) + exp(bbgr6.sim) )



# Irregular Reorganization Entry
Chyp   <- cbind(1,   0, log(tenure.range), 0, 0, 0 ,0,   mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) , 0*log(tenure.range),  0*log(tenure.range),   0*log(tenure.range), 0*log(tenure.range),   0*log(tenure.range))


cbgr2.sim <- Chyp%*% t(betadraw3[,1:23])
cbgr3.sim <- Chyp%*% t(betadraw3[,24:46])
cbgr4.sim <- Chyp%*% t(betadraw3[,47:69])
cbgr5.sim <- Chyp%*% t(betadraw3[,70:92])
cbgr6.sim <- Chyp%*% t(betadraw3[,93:115])


prob2cgr.sim <- exp(cbgr2.sim)/ (exp(0) + exp(cbgr2.sim) + exp(cbgr3.sim) + exp(cbgr4.sim) + exp(cbgr5.sim) + exp(cbgr6.sim) )
prob3cgr.sim <- exp(cbgr3.sim)/ (exp(0) + exp(cbgr2.sim) + exp(cbgr3.sim) + exp(cbgr4.sim) + exp(cbgr5.sim) + exp(cbgr6.sim) )
prob4cgr.sim <- exp(cbgr4.sim)/ (exp(0) + exp(cbgr2.sim) + exp(cbgr3.sim) + exp(cbgr4.sim) + exp(cbgr5.sim) + exp(cbgr6.sim) )
prob5cgr.sim <- exp(cbgr5.sim)/ (exp(0) + exp(cbgr2.sim) + exp(cbgr3.sim) + exp(cbgr4.sim) + exp(cbgr5.sim) + exp(cbgr6.sim) )
prob6cgr.sim <- exp(cbgr6.sim)/ (exp(0) + exp(cbgr2.sim) + exp(cbgr3.sim) + exp(cbgr4.sim) + exp(cbgr5.sim) + exp(cbgr6.sim) )




# First-Election Entry (regular)
Dhyp   <- cbind(1,   0, log(tenure.range), 0, 0, 1 ,0,   mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) , 0*log(tenure.range),  0*log(tenure.range),   0*log(tenure.range), 1*log(tenure.range),   0*log(tenure.range))

dbgr2.sim <- Dhyp%*% t(betadraw3[,1:23])
dbgr3.sim <- Dhyp%*% t(betadraw3[,24:46])
dbgr4.sim <- Dhyp%*% t(betadraw3[,47:69])
dbgr5.sim <- Dhyp%*% t(betadraw3[,70:92])
dbgr6.sim <- Dhyp%*% t(betadraw3[,93:115])


prob2dgr.sim <- exp(dbgr2.sim)/ (exp(0) + exp(dbgr2.sim) + exp(dbgr3.sim) + exp(dbgr4.sim) + exp(dbgr5.sim) + exp(dbgr6.sim) )
prob3dgr.sim <- exp(dbgr3.sim)/ (exp(0) + exp(dbgr2.sim) + exp(dbgr3.sim) + exp(dbgr4.sim) + exp(dbgr5.sim) + exp(dbgr6.sim) )
prob4dgr.sim <- exp(dbgr4.sim)/ (exp(0) + exp(dbgr2.sim) + exp(dbgr3.sim) + exp(dbgr4.sim) + exp(dbgr5.sim) + exp(dbgr6.sim) )
prob5dgr.sim <- exp(dbgr5.sim)/ (exp(0) + exp(dbgr2.sim) + exp(dbgr3.sim) + exp(dbgr4.sim) + exp(dbgr5.sim) + exp(dbgr6.sim) )
prob6dgr.sim <- exp(dbgr6.sim)/(exp(0) + exp(dbgr2.sim) + exp(dbgr3.sim) + exp(dbgr4.sim) + exp(dbgr5.sim) + exp(dbgr6.sim) )




# Regular Election Entry
Ehyp   <- cbind(1,   0, log(tenure.range), 0, 0, 0 ,1,   mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) , 0*log(tenure.range),  0*log(tenure.range),   0*log(tenure.range), 0*log(tenure.range),   1*log(tenure.range))

ebgr2.sim <- Ehyp%*% t(betadraw3[,1:23])
ebgr3.sim <- Ehyp%*% t(betadraw3[,24:46])
ebgr4.sim <- Ehyp%*% t(betadraw3[,47:69])
ebgr5.sim <- Ehyp%*% t(betadraw3[,70:92])
ebgr6.sim <- Ehyp%*% t(betadraw3[,93:115])


prob2egr.sim <- exp(ebgr2.sim)/ (exp(0) + exp(ebgr2.sim) + exp(ebgr3.sim) + exp(ebgr4.sim) + exp(ebgr5.sim) + exp(ebgr6.sim) )
prob3egr.sim <- exp(ebgr3.sim)/ (exp(0) + exp(ebgr2.sim) + exp(ebgr3.sim) + exp(ebgr4.sim) + exp(ebgr5.sim) + exp(ebgr6.sim) )
prob4egr.sim <- exp(ebgr4.sim)/ (exp(0) + exp(ebgr2.sim) + exp(ebgr3.sim) + exp(ebgr4.sim) + exp(ebgr5.sim) + exp(ebgr6.sim) )
prob5egr.sim <- exp(ebgr5.sim)/ (exp(0) + exp(ebgr2.sim) + exp(ebgr3.sim) + exp(ebgr4.sim) + exp(ebgr5.sim) + exp(ebgr6.sim) )
prob6egr.sim <- exp(ebgr6.sim)/(exp(0) + exp(ebgr2.sim) + exp(ebgr3.sim) + exp(ebgr4.sim) + exp(ebgr5.sim) + exp(ebgr6.sim) )




## Coalition-Competing Exit (with Irregular Reorganization Entry as base category)

# Pr(Coalition Competing|Rebel entry) - Pr(Coalition Competing|Irregular Reorganization Entry)

Diff.P2.sim.ac <-  prob2agr.sim - prob2cgr.sim

# Pr(Coalition Competing|Regime-Change entry) - Pr(Coalition Competing|Reorganization Entry)

Diff.P2.sim.bc <-  prob2bgr.sim - prob2cgr.sim




# Pr(Coalition Competing|First Regular election entry) - Pr(Coalition Competing|Reorganization entry)

Diff.P2.sim.dc <- prob2dgr.sim - prob2cgr.sim

# Pr(Coalition Competing|Regular Election entry) - Pr(Coalition Competing|Reorganization eentry)

Diff.P2.sim.ec <- prob2egr.sim - prob2cgr.sim



##
## Figure 3
##

## Effect of Rebel & Regime-Change Coup Entries on Coalition-Competing Exit -- Regular Entry as the base category)


windows()
par(mar=c(4,4,3,1), mfrow=c(2   ,2))
# A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot.
# The default is c(5, 4, 4, 2) + 0.1.


plot(tenure.range, rowMeans(Diff.P2.sim.gz) ,type="l", xlim=c(1, 10),  xaxt='n', ylim=c(-0.02, 0.03), main="Effect of Rebel Entry (v.s. Regular)", lwd=3,ylab="Pr(Competing Exit|Rebel Entry) - Pr(Competing Exit|Regular Entry)", xlab="Tenure")
lines(tenure.range, t(apply(Diff.P2.sim.gz,1,sort))[,9750],type="l",lty=3,lwd=2    )
lines(tenure.range,  t(apply(Diff.P2.sim.gz,1,sort))[,250],type="l",lty=3,lwd=2  )
abline(h=0,lty=3 )
axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T,  line=-2)
legend(1, 0.03,  bty="n",   c( "Effect of Rebel Entry (v.s. Regular)" , "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")



plot(tenure.range, rowMeans(Diff.P2.sim.yz) ,type="l", xlim=c(1, 10), xaxt='n', ylim=c(-0.02 , 0.03), main="Effect of Regime-Change Coup Entry (v.s. Regular)",   lwd=3,ylab="Pr(Competing Exit|Reg-Change Coup) - Pr(Competing Exit|Regular Entry)", xlab="Tenure")
lines(tenure.range, t(apply(Diff.P2.sim.yz,1,sort))[,9750],type="l",lty=3,lwd=2   )
lines(tenure.range,  t(apply(Diff.P2.sim.yz,1,sort))[,250],type="l",lty=3,lwd=2  )
#legend(x=500, y=0.6, bty="n", legend="Incumbency Advantage", lty=1, lwd=3, seg.len=0.75, col="blue", x.intersp=0.5)
abline(h=0,lty=3 )

axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)
legend(1, 0.03,  bty="n",   c( "Effect of Regime-Change Coup Entry (v.s. Regular)" , "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")



plot(tenure.range, rowMeans(Diff.P2.sim.ac ) ,type="l", xlim=c(1, 10),  xaxt='n', ylim=c(-0.035, 0.05),  main="Effect of Rebel Entry (v.s. Reorganization)", lwd=3,ylab="Pr(Competing Exit|Rebel Entry) - Pr(Competing Exit|Reorganization Entry)", xlab="Tenure")
lines(tenure.range, t(apply(Diff.P2.sim.ac ,1,sort))[,9750],type="l",lty=3,lwd=2    )
lines(tenure.range,  t(apply(Diff.P2.sim.ac ,1,sort))[,250],type="l",lty=3,lwd=2  )
abline(h=0,lty=3 )
axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T,  line=-2)
legend(1, 0.05,  bty="n",   c( "Effect of Rebel Entry (v.s. Irregular Reorganization)" , "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")



plot(tenure.range, rowMeans(Diff.P2.sim.bc ) ,type="l", xlim=c(1, 10), xaxt='n', ylim=c(-0.035, 0.05), main="Effect of Regime-Change Coup Entry (v.s. Reorganization)",    lwd=3,ylab="Pr(Competing Exit|Reg-Change Coup) - Pr(Competing Exit|Reorganization Entry)", xlab="Tenure")
lines(tenure.range, t(apply(Diff.P2.sim.bc ,1,sort))[,9750],type="l",lty=3,lwd=2   )
lines(tenure.range,  t(apply(Diff.P2.sim.bc  ,1,sort))[,250],type="l",lty=3,lwd=2  )
#legend(x=500, y=0.6, bty="n", legend="Incumbency Advantage", lty=1, lwd=3, seg.len=0.75, col="blue", x.intersp=0.5)
abline(h=0,lty=3 )

axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)
legend(1, 0.05,  bty="n",   c( "Effect of Regime-Change Coup Entry (v.s. Irregular Reorganization)" , "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")



# save

savePlot(filename="C:/Users/pjb13215/Dropbox/Entry & Exit manner in Dictatorships/ISQ Submission/R&R/dis_meffect_competing_removal", type = "pdf")






##
## Figures in Appendix (Disaggregated Reorganization Entry)

windows()
par(mar=c(4,4,3,1), mfrow=c(2,2))
# A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot.
# The default is c(5, 4, 4, 2) + 0.1.

plot(tenure.range, rowMeans(Diff.P2.sim.rz) ,type="l", xlim=c(1, 10), xaxt='n', ylim=c(-0.01, 0.08), lwd=3,ylab="Pr(Competing Exit|Reshuffling-Coup) - Pr(Competing Exit|Regular)", xlab="Tenure")
lines(tenure.range, t(apply(Diff.P2.sim.rz,1,sort))[,9750],type="l",lty=3,lwd=2    )
lines(tenure.range,  t(apply(Diff.P2.sim.rz,1,sort))[,250],type="l",lty=3,lwd=2  )
abline(h=0,lty=3 )
axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T,  line=-2)
legend(1, 0.08,  bty="n",   c( "Effect of Reshuffling-Coup Entry (v.s. Regular)" , "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")

plot(tenure.range, rowMeans(Diff.P2.sim.oz) ,type="l", xlim=c(1, 10) , xaxt='n', ylim=c(-0.02 , 0.08),   lwd=3,ylab="Pr(Competing Exit|Forced Selection) - Pr(Competing Exit|Regular)", xlab="Tenure")
lines(tenure.range, t(apply( Diff.P2.sim.oz,1,sort))[,9750],type="l",lty=3,lwd=2   )
lines(tenure.range,  t(apply(Diff.P2.sim.oz,1,sort))[,250],type="l",lty=3,lwd=2  )
#legend(x=500, y=0.6, bty="n", legend="Incumbency Advantage", lty=1, lwd=3, seg.len=0.75, col="blue", x.intersp=0.5)
abline(h=0,lty=3 )

axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)
legend(1, 0.08,  bty="n",   c( "Effect of Forced Selection Entry (v.s. Regular)" , "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")



plot(tenure.range, rowMeans(Diff.P2.sim.dc) ,type="l", xlim=c(1, 10), xaxt='n', ylim=c(-0.03, 0.04), lwd=3,ylab="Pr(Competing Exit|First-Election) - Pr(Competing Exit|Reorganization)", xlab="Tenure")
lines(tenure.range, t(apply(Diff.P2.sim.dc,1,sort))[,9750],type="l",lty=3,lwd=2    )
lines(tenure.range,  t(apply(Diff.P2.sim.dc,1,sort))[,250],type="l",lty=3,lwd=2  )
abline(h=0,lty=3 )
axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T,  line=-2)
legend(1, 0.04,  bty="n",   c( "Effect of First-Election Entry (v.s. Irregular Reorganization)" , "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")

plot(tenure.range, rowMeans(Diff.P2.sim.ec) ,type="l", xlim=c(1, 10) , xaxt='n', ylim=c(-0.03 , 0.04),   lwd=3,ylab="Pr(Competing Exit|Regular Election) - Pr(Competing Exit|Reorganization)", xlab="Tenure")
lines(tenure.range, t(apply( Diff.P2.sim.ec,1,sort))[,9750],type="l",lty=3,lwd=2   )
lines(tenure.range,  t(apply(Diff.P2.sim.ec,1,sort))[,250],type="l",lty=3,lwd=2  )
#legend(x=500, y=0.6, bty="n", legend="Incumbency Advantage", lty=1, lwd=3, seg.len=0.75, col="blue", x.intersp=0.5)
abline(h=0,lty=3 )

axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)
legend(1, 0.04,  bty="n",   c( "Effect of Regular Election Entry (v.s. Irregular Reorganization)" , "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")






# save
savePlot(filename="C:/Users/pjb13215/Dropbox/Entry & Exit manner in Dictatorships/ISQ Submission/R&R/dis_meffect_competing_removal_app", type = "pdf")










##
## Figures 4 & 5 (Logit Model )
##

#DV: Regime- coup exit

betadraw <- read.table("Data/VCVmatrix_v3.csv", header=TRUE, sep=",")  # this is a 10000 draws of simulated coefficients.
dim(betadraw) # 10000 * 12



f <-  regimecoup_exit ~ regimecoup_entry*log(tenure) + leadercoup_entry*log(tenure) + lngdp + growth + military + monarch + party + log_milper +
 age  +  dwin + dlose + ddraw + incidencev414

model1 <- glm(f, family=binomial(link = "logit"), data=dat)

coefficients(model1)
tenure.range <- seq(1, 10, .1)


# Regime-Replacement Coup entry leader

Z <- cbind(1,   1, log(tenure.range), 0,    mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 1*log(tenure.range),  0*log(tenure.range))

coefficients(model1)
Zb <- Z%*%coefficients(model1)
ZP <- exp(Zb)/(1+exp(Zb))      # 201*1 vector of P


Zb.sim <- Z %*% t(betadraw)
ZP.sim <- exp(Zb.sim)/(1+exp(Zb.sim))   # 201*1000 simulated probabilities.
ZP.sim.graph <- t(apply(ZP.sim,1,sort)) # 201*1000 simulated probabilities, sorted to identify the confidence bounds.



# Leader Reshuffling Coup entry leader

X <- cbind(1,   0, log(tenure.range), 1,    mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE),
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 0*log(tenure.range),  1*log(tenure.range))

coefficients(model1)
Xb <- X%*%coefficients(model1)
XP <- exp(Xb)/(1+exp(Xb))      # 201*1 vector of P

Xb.sim <- X %*% t(betadraw)
XP.sim <- exp(Xb.sim)/(1+exp(Xb.sim))   # 201*1000 simulated probabilities.
XP.sim.graph <- t(apply(XP.sim,1,sort)) # 201*1000 simulated probabilities, sorted to identify the confidence bounds.



# No coup-entry leader

W <- cbind(1,   0, log(tenure.range),0,    mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE),
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 0*log(tenure.range),  0*log(tenure.range))


Wb <- W %*%coefficients(model1)
WP <- exp(Wb)/(1+exp(Wb))      # 201*1 vector of P

Wb.sim <- W %*% t(betadraw)
WP.sim <- exp(Wb.sim)/(1+exp(Wb.sim))   # 201*1000 simulated probabilities.
WP.sim.graph <- t(apply(WP.sim,1,sort)) # 201*1000 simulated probabilities, sorted to identify the confidence bounds.


Diff.P.sim.zw <- ZP.sim - WP.sim
Diff.P.sim.xw <- XP.sim - WP.sim


#
# Figure 4
#

windows()
par(mar=c(5,5,4,3))
# A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot.
# The default is c(5, 4, 4, 2) + 0.1.
par(mfrow=c(1,2))


plot(tenure.range, rowMeans(Diff.P.sim.zw) ,type="l", ylim=c(-0.01, 0.015) ,xlim=c(1, 10), lwd=3,xaxt='n', ylab="Pr(Reg-Change Exit|Reg-Change Entry)- Pr(Reg-Change Exit|Non-Coup Entry)", xlab="Tenure")
abline(h=0)
lines(tenure.range, t(apply(Diff.P.sim.zw,1,sort))[,9750],type="l",lty=3,lwd=2)
lines(tenure.range,  t(apply(Diff.P.sim.zw,1,sort))[,250],type="l",lty=3,lwd=2)
abline(h=0,lty=3 )
axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")


legend(1, 0.015,   bty="n",  c("Effect of Regime-Change Coup Entry (v.s. Non-Coup Entry)", "95 % CI"),lwd=c(3,2), lty=c(1,3) , col="black")
#mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)



plot(tenure.range, rowMeans(Diff.P.sim.xw) ,type="l", ylim=c(-0.01, 0.02) ,xlim=c(1, 10), lwd=3,xaxt='n', ylab="Pr(Reg-Change Exit|Reshuffling Entry)- Pr(Reg-Change Exit|Non-Coup Entry)", xlab="Tenure")
abline(h=0)
lines(tenure.range, t(apply(Diff.P.sim.xw,1,sort))[,9750],type="l",lty=3,lwd=2)
lines(tenure.range,  t(apply(Diff.P.sim.xw,1,sort))[,250],type="l",lty=3,lwd=2)
abline(h=0,lty=3 )
axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")


legend(1, 0.02, bty="n",  c("Effect of Reshuffling Coup Entry (v.s. Non-Coup Entry)", "95 % CI"),lwd=c(3,2), lty=c(1,3) , col="black")
#mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)




# save

savePlot(filename="C:/Users/pjb13215/Dropbox/Entry & Exit manner in Dictatorships/ISQ Submission/R&R/meffect_regime_coup", type = "pdf")









##
## DV: Leader Reshuffle coup exit

betadraw <- read.table("Data/VCVmatrix_v4.csv", header=TRUE, sep=",")  # this is a 10000 draws of simulated coefficients.
dim(betadraw) # 10000 * 12


h <-  leadercoup_exit ~ regimecoup_entry*log(tenure) + leadercoup_entry*log(tenure) + lngdp + growth + military + monarch + party + log_milper +
 age  +  dwin + dlose + ddraw + incidencev414

model2 <- glm(h, family=binomial(link = "logit"), data=dat)

tenure.range <- seq(1, 10, .1)


# Regime-Replacement Coup entry leader

Z <- cbind(1,   1, log(tenure.range), 0,    mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 1*log(tenure.range),  0*log(tenure.range))

coefficients(model2)
Zlb <- Z%*%coefficients(model2)
ZlP <- exp(Zlb)/(1+exp(Zlb))      # 201*1 vector of P

Zlb.sim <- Z %*% t(betadraw)
ZlP.sim <- exp(Zlb.sim)/(1+exp(Zlb.sim))   # 201*1000 simulated probabilities.
ZlP.sim.graph <- t(apply(ZlP.sim,1,sort)) # 201*1000 simulated probabilities, sorted to identify the confidence bounds.


# Leader Reshuffling Coup entry leader

X <- cbind(1,   0, log(tenure.range), 1,    mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 0*log(tenure.range),  1*log(tenure.range))

Xlb <- X%*%coefficients(model2)
XlP <- exp(Xlb)/(1+exp(Xlb))      # 201*1 vector of P

Xlb.sim <- X %*% t(betadraw)
XlP.sim <- exp(Xlb.sim)/(1+exp(Xlb.sim))   # 201*1000 simulated probabilities.
XlP.sim.graph <- t(apply(XlP.sim,1,sort)) # 201*1000 simulated probabilities, sorted to identify the confidence bounds.


# No coup-entry leader

W <- cbind(1,   0, log(tenure.range),0,    mean( lngdp , na.rm=TRUE),  mean(growth, na.rm=TRUE), median(military, na.rm=TRUE) ,median(monarch, na.rm=TRUE) , median(party, na.rm=TRUE) ,
 mean(log_milper , na.rm=TRUE),  mean(age , na.rm=TRUE),  mean(dwin, na.rm=TRUE) ,mean( dlose , na.rm=TRUE) , mean( ddraw, na.rm=TRUE) , median( incidencev414, na.rm=TRUE) ,
 0*log(tenure.range),  0*log(tenure.range))


Wlb <- W %*%coefficients(model2)
WlP <- exp(Wlb)/(1+exp(Wlb))      # 201*1 vector of P

Wlb.sim <- W %*% t(betadraw)
WlP.sim <- exp(Wlb.sim)/(1+exp(Wlb.sim))   # 201*1000 simulated probabilities.
WlP.sim.graph <- t(apply(WlP.sim,1,sort)) # 201*1000 simulated probabilities, sorted to identify the confidence bounds.





Diff.lP.sim.zw <- ZlP.sim - WlP.sim
Diff.lP.sim.xw <- XlP.sim - WlP.sim


#
# Figure 5
#
windows()
par(mar=c(5,5,4,3))
# A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot.
# The default is c(5, 4, 4, 2) + 0.1.
par(mfrow=c(1,2))


plot(tenure.range, rowMeans(Diff.lP.sim.zw ) ,type="l", ylim=c(-0.01, 0.05) ,xlim=c(1, 10), lwd=3,xaxt='n', ylab="Pr(Reg-Change Exit|Reg-Change Entry)- Pr(Reg-Change Exit|Non-Coup Entry)", xlab="Tenure")
abline(h=0)
lines(tenure.range, t(apply(Diff.lP.sim.zw ,1,sort))[,9750],type="l",lty=3,lwd=2)
lines(tenure.range,  t(apply(Diff.lP.sim.zw ,1,sort))[,250],type="l",lty=3,lwd=2)
abline(h=0,lty=3 )
axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")


legend(1, 0.05,   bty="n",  c("Effect of Regime-Change Coup Entry (v.s. Non-Coup Entry)", "95 % CI"),lwd=c(3,2), lty=c(1,3) , col="black")
#mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)



plot(tenure.range, rowMeans(Diff.lP.sim.xw ) ,type="l", ylim=c(-0.01, 0.06) ,xlim=c(1, 10), lwd=3,xaxt='n', ylab="Pr(Reg-Change Exit|Reshuffling Entry)- Pr(Reg-Change Exit|Non Coup Entry)", xlab="Tenure")
abline(h=0)
lines(tenure.range, t(apply(Diff.lP.sim.xw ,1,sort))[,9750],type="l",lty=3,lwd=2)
lines(tenure.range,  t(apply(Diff.lP.sim.xw ,1,sort))[,250],type="l",lty=3,lwd=2)
abline(h=0,lty=3 )
axis(1, at=c(1, 5, 10), label=c(1, 5, 10), col.axis="black")


legend(1, 0.06, bty="n", c("Effect of Reshuffling Coup Entry (v.s. Non-Coup Entry)", "95 % CI"),lwd=c(3,2), lty=c(1,3) , col="black")
#mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)




# save

savePlot(filename="C:/Users/pjb13215/Dropbox/Entry & Exit manner in Dictatorships/ISQ Submission/R&R/meffect_reshuffling_coup", type = "pdf")














# end of file.
